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Abstract: 

The  combination  of  data  and  information  from  multiple  sources  or  sensors 
constitutes  one  of  the  most  important  problems  in  signal  processing  to-day.  This  process 
of  data  fusion  has  been  defined  as  one  of  dealing  with  the  association,  correlation  and 
combination  of  data.  In  this  report,  an  exact  maximum  likelihood  method  is  used  to  solve 
the  problem  of  registration  errors  in  data  fusion.  Registration  is  defined  as  the  co-ordinate 
conversion  of  multiple  source  data.  Error-free  registration  is  a  pre-requisite  in  the  process 
of  data  fusion.  In  the  exact  maximum  likelihood  method  (EML),  a  likelihood  function  is 
constructed  based  on  the  errors  of  co-ordinate  transformation  of  the  local  sensor  locations 
to  a  common  system  plane.  Optimization  is  then  carried  out  in  a  recursive  two-step  Gauss- 
Newton  type  procedure,  which  produces  the  correct  system  biases  for  proper  registration. 
Simulation  results  show  that  EML  is  more  efficient  when  compared  with  conventional 
methods  of  registration. 


P5G0711.PDF  [Page:  7  of  48] 


t 


Introduction 

Data  fusion  has  been  defined  as  "a  process  dealing  with  association,  correlation 
and  combination  of  data  and  information  from  single  and  multiple  sources  to  achieve 
refined  position  and  identity  estimates"  [1].  This  aspect  of  signal  processing  has  increased 
in  importance  in  recent  years  [2]  owing  to  the  advance  in  technology  of  sensor  systems 
which  are  becoming  more  sophisticated  and  complex.  As  more  of  these  sensor  systems 
become  available,  the  combination  of  information  from  a  multitude  of  sources  and  the 
decision  on  proper  reaction  turn  into  a  significant  problem  [3].  This  fusion  process  is 
usually  implemented  at  three  levels.  Level  1  includes  registration,  association,  filtering  and 
identification.  Level  2  refers  to  situation  assessment  and  level  3  determines  whether  a 
threat  is  present  [4].  In  this  report,  the  focus  is  on  level  1  processing  and,  in  particular,  on 
the  registration  problem  [5]. 

The  process  of  proper  co-ordinate  conversion  of  multiple  source  data  is  called 
registration.  Generally,  this  is  required  because  the  sensors  are  not  aligned  properly  for 
fusion.  The  major  source  of  registration  error  are  the  azimuth  bias  error  with  respect  to  a 
common  reference  and  the  range  offset  errors  of  each  sensor.  These  errors  in  turn  will 
introduce  biases  in  fusion  and  may  generate  "ghost"  targets  in  multisensor  signal 
processing  [6,7].  To  "register"  a  sensor  system,  each  sensor  has  to  be  initialized  and  then 
aligned  to  a  common  reference  frame.  Initialization  refers  to  the  independent  localization 
of  each  sensor  in  the  system  co-ordinates.  Relative  alignment  is  done  by  using  common 
targets  and  collecting  data  points  from  each  sensor.  A  set  of  system  errors  (for  proper 
registration)  is  then  computed  and  used  to  adjust  subsequent  incoming  data  in  the  data 
fusion  process. 

Several  registration  methods  have  been  proposed.  These  include  the  simple 
averaging  method  called  Real  Time  Quality  Control  (RTQC)  [8],  the  least  squares  (LS) 
method  [9],  the  maximum  likelihood  (ML)  solution  [10]  and  the  generalized  least  squares 
(GLS)  technique  [5].  The  RTQC  routine  computes  the  system  errors  by  averaging  the 
measurement  data.  The  LS  formulates  the  registration  as  an  ordinary  or  unweighted  least 
squares  problem,  while  GLS  includes  the  target  and  observation  error  covariance  in  the 
formulation.  These  methods  are  found  to  yield  reasonable  results  only  when  the 
measurement  noise  is  relatively  small.  The  major  drawback  is  the  assumption  that  the 
measurements  from  different  sensors  are  "equal"  in  a  common  co-ordinate  system.  Thus, 
the  measurement  noise  from  individual  sensors  is  not  taken  into  account  properly. 


This  report  proposes  an  exact  maximum  likelihood  method  for  registration.  The 
likelihood  function  is  derived  directly  from  the  measurement  data,  which  are  assumed  to 
contain  two  sets  of  parameters:  the  true  target  positions  and  the  system  errors.  These  two 
sets  of  variables  are  found  to  be  partially  separable  in  the  likelihood  function,  and  thus  a 
two-step  recursive  optimization  procedure  is  possible.  An  approximation  to  the  Hessian  is 
also  made  to  ensure  convergence  of  the  optimization  procedure.  In  the  following  sections, 


P5Q0711.PDF  [Page:  8  of  48] 


5 


we  present  the  conventional  formulation  of  the  registration  problem,  the  exact  maximum 
likelihood  formulation,  the  algorithm,  computer  simulation  and  real  data  analysis.  Some 
conclusions  are  then  drawn. 


The  Registration  Problem 

Consider  two  sensors  A  and  B  which  measure  the  range  and  azimuth  of  a  common 
target.  Without  loss  of  generality,  let  A  be  located  at  the  origin  and  B  at  (u,v)  as  shown  in 
Figure  1.  Let  Tk  denote  the  k-th  target,  and  let  {rA(k),  0A(k)>,{rB(k),  0B(k)>  represent 
the  range  and  azimuth  of  the  k-th  target  as  measured  by  A  and  B  respectively.  The 
corresponding  sensor  biases  (system  errors)  are  denoted  by  {ArA,  A0a>  and  {ArB,  A0B>. 
If  {x'(k),  y'(k)>  is  the  true  position  for  the  target  Tk ,  we  have, 


x'(k)  -  [rA(k)  -  ArA]  sin[0A(k)  -  A0a] 

y'(k)  -  [rA(k)  -  ArA]  cos[0 A(k)  -  A0 A  ] ,  ( 1 ) 

for  sensor  A,  and 

x'(k)  -  [rB(k)  -  ArB]  sin[0B(k)  -  A0B]  +  u 

y'(k)  =  [rB(k)  -  ArB]  cos[0B(k)  -  A0B  ]  +  v  ,  (2) 

for  sensor  B.  These  two  sets  of  equations  can  be  combined  and  be  represented  in  matrix 
form,  giving  us, 

L(k)  H  =  Ax  (k)  (3) 

where  u  -  [  A0a,  ArA,  A0B,  ArB]T  , 


L(k)  «= 


^rA(k)cos0A(k)  sin0A(k) 
^-rA(k)sin0A(k)  cos0A(k) 


-rB(k)cos0B(k) 

rB(k)sin0B(k) 


-sin0B(k)  N 
-cos0B(k), 


and 


Ax  (k)  - 


f  rA(k)sin0A(k)-rB(k)sin0B(k)-u  N 
^rA(k)cos0  A(k)-rB(k)cos0B(k)-v  , 


(4) 

(5) 


Equation  (3)  is  the  basic  registration  equation  for  all  conventional  methods.  The 
only  difference  is  the  way  n  is  solved.  The  RTQC  routine  separates  the  measurements 
into  two  groups  by  the  site  line  joining  the  sensors  A  and  B.  It  then  averages  the 
measurements  in  each  group,  and  the  registration  error  vector  3}  can  then  be  obtained  by 
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solving  a  4x4  matrix  [8].  The  LS  method  solves  for  U  by  finding  the  least  squares  solution 
based  on  K  measurements.  Thus,  equation  (3)  becomes, 

L  H  -  As  (6) 

where  L  -  [L(1)T,  L(2)T, . L(K)T]T ,  and  As  -  [A*T(1),  AsT(2), ...,  A*T(K)]T ,  and 

Dls  ”  (LTL)'1LT  Ax .  (7) 

The  LS  solution  is  based  on  the  assumption  that  all  measurements  are  equally 
weighted  [10].  The  ML  and  the  GLS  methods  are  introduced  to  incorporate  the 
covariances  of  the  measurement  noise  to  weight  the  measurements.  Since  these  two 
methods  do  not  seem  to  have  much  significant  improvement  over  the  LS  technique, 
especially  in  the  real-life  applications  [11],  our  comparisons  will  be  made  with  the  LS 
method  only. 

The  Exact  Maximum  Likelihood  Method 

The  limitations  of  the  conventional  formulation  for  the  registration  errors  as  shown 
by  equation  (3)  are  evident  when  measurement  noise  are  present  in  the  sensors.  The 
formulation,  whether  solved  by  simple  averaging  or  least  squares,  assumes  no  noise  in  the 
measured  data,  as  shown  in  equations  (1)  and  (2).  Any  subsequent  correction  using  the 
computed  biases  in  range  and  azimuth  will  further  enhance  the  random  error  in  the 
measurements.  This  explains  the  limited  success  conventional  techniques  had  with 
registration. 

In  the  EML  method,  the  random  measurement  noise  of  the  sensors  are  included  in 
the  formulation.  Let  {xA(k),  yA(k)>  and  <xB(k),  yB(k)>  denote  the  Cartesian  co-ordinates 
of  the  target  Tk  in  the  system  plane  as  measured  by  sensors  A  and  B  respectively.  The 
actual  range  and  azimuth  measured  by  A  and  B  are  {r'A(k),  0'A(k)  >  and  {r'B(k),  0'B(k)}. 
Then  using  nj(k)  as  the  random  noise  components,  we  have 

xA(k)  -  [r'A(k)  +  ArA]  sin[0'A(k)+AeA]  +  n^k) 

y A00  “  [r’A(k)  +  ArA]  cos[0'A(k)+A0A]  +  n2(k),  (8) 

and 

xB(k)  “  [r'B(k)  +  ArB]  sin[0'B(k)+A0B]  +  u  +  n3(k) 

yB(k)  -  [r'B(k)  +  ArB]  cos[0'B(k)+A0B]  +  v  +  n4(k).  (9) 
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Generally,  the  system  biases  A0a,  Ata,  A0b,  ArB  are  small  compared  to  measured 
values.  The  equations  can  be  well  approximated  by  retaining  first  order  terms  in  the 
biases.  Thus, 

xA(k)  =  [r'A(k)  +  ArA]sin0'A(k)  +  r'A(k)cos0'A(k)A0A  +  n^k), 
yA(k)  =  [r'A(k)  +  ArA]cos0'A(k)  -  r'A(k)sin0,A(k)A0A  +  n2(k), 
xB(k)  =  [r'B(k)  +  ArB]sin0'B(k)  +  r'B(k)cos0'B(k)A0B  +  u  +  n3(k), 
yB(k)  =  [r'A(k)  +  ArA]cos0'B(k)  -  r'B(k)sin0'B(k)A0B  +  v  +  n4(k).  (10) 

Now  when  {x'(k),  y'(k)}  is  the  actual  position  of  Tj.  in  the  system  plane,  we  have 
x'(k)  -  r'A(k)  sin0'A(k)  -  r'B(k)  sin0'B(k)  +  u, 

y'(k)  »  r'A(k)  cos0'A(k)  »  r'B(k)  cos0'B(k)  +  v.  (1 1) 

These  can  be  substituted  into  (10)  giving, 

xA(k)  =  x'(k)  +  x'(k)  +  A0Ay'(k)  +  n^k), 

yA(k)  -  y'(k)  +  y  W  -  A0Ax'(k)  +  n2(k), 

xB(k)  -  x'(k)  +  x'(k)  +  A0By'(k)  -  u  -  v  A0B  +  n3(k), 

yB(k)  -  y'(k)  +  y'(k)  +  A0Bx'(k)  -  v  +  u  A0B  +  n4(k).  (12) 

In  matrix  form,  equation  (12)  becomes, 
x(k)  =  A(k)  2}  +  b(k)  +  n(k) 


where 


x(k)  =  [xA(k),  yA(k),  xB(k),  yB(k)]T, 
b(k)  =  [x'(k),  y'(k),  x'(k),  y'(k)]T, 


(13) 
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n(k)  -  [ni(k),  n2(k),  n3(k),  n4(k)]T, 


H  -  [A9a,  ArA,  A0b,  ArB]T, 


(14) 


and  the  matrix  A(k)  is  a  block  diagonal  matrix  given  by  A(k)  -  diag[An(k),  A22(k)], 
where  the  blocks  are  defined  as 


Au(k) 


^  y'(k) 
-x'(k) 


x'(k)  A 

r’A(k) 

m 

OOj 


and 


A22(k)  — 


^  yW-v 
-x'(k)+u 


x'(k)-uA 

1*800 

y'.(£hx 


(15) 


2 

n(k)  is  the  random  measurement  noise  vector  with  a  covariance  matrix  given  by  on  I.  i(k), 
rj  are  the  target  measurement  vector  and  the  system  bias  vector  for  the  target  Tk  .  The 
matrix  A(k)  and  the  vector  Mk)  are  independent  of  the  system  biases  and  are  determined 
by  the  actual  location  of  the  target  in  the  system  plane. 

Assuming  {n(k)}  to  be  independently  and  normally  distributed,  the  conditional 
probability  density  function  of  the  measurements (x(k),  k-l,2,....,K>  is, 


p([x(l),i(2),x(3),....^c(K)]) 


=nk_i — l—  exp  [i(k)  -  A(k)n  -  Mk)]T[2c(k)  -  A(k)u  -  h(k)]  >  (16) 

(2n)2on  2  °n 

The  corresponding  negative  log  likelihood  function,  without  the  constant  factor  has  the 
form. 


K 

J  -  -  log  p  -  2K  log  (2nol)  +  S II  x(k)  -  A(k)u  -  b(k)llF2.  (17) 

2on  k-1 

where  II  •  llF  denotes  the  Frobenius  norm.  Based  on  the  principle  of  maximum  likelihood 
estimation,  the  estimates  of  the  unknown  parameters  are  obtained  by  the  minimizing 

arguments  of  (17)  [12].  Let  £(k)  -  [x'(k),  y'(k)]T  be  the  actual  position  vector  of  the 

2  . 

target  in  the  system  plane.  Fixing  ^(k)  and  n ,  while  minimizing  J  with  respect  to  an  gives 
the  estimate  for  the  noise  variance: 

£  II  x(k)  -  A(k)n  -  £>(k)llF2. 

k-l 


(18) 
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Substituting  this  back  into  (17),  the  estimates  for  £(k)  and  u  are  obtained  as  the 
minimizing  arguments  for  the  function  J: 

[£(k),  n  ]  -  arg  min  J,  (19) 

where 

J-4  him- A(k)n  -  b(k)llF2.  (20) 

k-l 


In  general,  J  is  a  nonlinear  function  of  the  two  parameter  vectors.  However,  it  is 
seen  that  the  parameters  £(k)  and  n  are  separable  in  (20),  and  an  alternating  optimization 
technique  can  be  applied  [13].  This  technique  iterates  between  two  steps,  each  minimizing 
one  of  the  parameter  vectors  until  convergence  is  reached. 

The  Optimization  Procedure 

Estimation  of  the  system  bias  vector 

The  bias  vector  r\  is  first  estimated  while  the  true  target  position  vector  ^(k)  is 
fixed.  This  means  that  the  estimate  is  the  solution  of, 


3J 

an 


--2 


2  AT(k)[x(k)  -  A(k)u  -  b(k)]  —  0. 


k-l 


(21) 


Solving  (21),  the  estimate  of  the  bias  vector  is  given  by, 

21-  <  f  AT(k)A(k)}-l  XAT(k)[x(k)-32(k)]. 
k-l  k-l 


(22) 


Estimating  the  true  target  position  vector 

Using  the  estimate  in  (22)  the  true  target  position  vector  can  be  calculated  as  the 
minimizing  arguments  of  the  likelihood  function  J.  It  can  be  seen  that  since  J  is  a  sum  over 
non-negative  terms,  each  of  which  relates  to  a  specific  k,  the  minimization  can  be 
performed  term-by-term  so  that  we  obtain, 

£(k)  ■=  arg  min  Jk  -=  arg  min  II  x(k)  -  A(k)ii  -  h(k)llF2.  (23) 

The  non-linear  problem  in  (23)  can  be  solved  numerically  by  using  a  Gauss-Newton  type 
method  [13].  The  method  can  be  represented  by  the  iterative  equation, 
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|(P+l)(k)-|(P)(k)-tlpHk‘Gk  (24) 

where  lip  is  the  p-th  iteration  step  length,  Hk  is  the  Hessian  matrix  for  the  function  Jk  with 
respect  to  the  parameters  to  be  estimated,  and  Gk  is  the  gradient  given  by, 

Gk  -  2  Rk  [x(k)  -  b(k)]  -  2Rky(k).  (25) 

Rk  is  the  Jacobian  matrix  of  ^(k)  with  respect  to  the  target  position  vectors  £(k)  and  can 
be  obtained  as, 


The  Hessian  of  Jk  can  be  computed  as  its  second  order  derivative  with  respect  to  the 
parameters  to  be  estimated.  Its  role  is  to  modify  the  gradient  in  order  to  achieve  faster 
convergence.  At  each  iteration,  the  Hessian  must  be  positive  definite  to  guarantee 
convergence.  In  practice,  the  true  second  order  derivative  is  not  guaranteed  to  be  positive 
definite  at  each  iteration.  To  avoid  non-convergence,  the  Hessian  is  approximated  by 

Hk  -  2  RkRkT  .  (27) 

Thus  the  Hessian  will  be  a  positive  semi-definite  matrix.  This  approximation  is  justified  by 
observing  that  the  system  errors  are  generally  much  smaller  than  the  range  measurements 
of  the  target. 

The  EML  registration  algorithm  can  be  summarized  as  follows. 


•  Set  a  threshold  e  and  set  the  iteration  index  p  to  zero. 

•  Obtain  an  initial  estimate  of  the  true  target  position  vectors  £(P>(k)  for  k  =  1,2,..., 

K,  and  compute  the  estimate  according  to  (22). 
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•  Minimize  the  function  Jk ,  using  the  approximate  Hessian,  to  obtain  the  (p+l)-th 
estimate  £(P+1)(k). 

•  Estimate  the  (p+l)-th  system  bias  vector  xj(P+1)  using  (22).  Compute  the 
difference  norm,  d  -  lli|(P+1)  -  ^(p)II2  and  check  it  against  the  threshold  e.  Repeat 
from  step  3  if  d  >  £  and  set  index  p  to  p+1;  otherwise  stop. 

We  note  that  the  last  two  steps  can  be  interchanged  in  the  sense  that  the  threshold 
can  be  applied  either  to  the  system  bias  vector  or  the  true  target  position  vector.  With  the 
Gauss-Newton  method,  where  the  approximated  Hessian  is  used,  each  separate  estimation 
of  the  system  bias  vector  and  of  the  target  position  vector  will  not  increase  the  objective 
function  Jk  .  Thus,  although  convergence  to  a  global  minimum  is  not  guaranteed  the 
algorithm  is  seen  to  be  convergent. 

The  algorithm  can  be  further  simplified,  since  the  matrix  A(k)  is  block  diagonal. 
This  means  the  size  of  the  matrices  to  be  inverted  is  reduced  by  a  factor  of  1/2. 


Simulation  and  real  data  analysis 

In  the  simulation  studies,  we  examine  the  performance  and  robustness  of  the  EML 
method  using  track  patterns  shown  in  Figure  2.  Sensor  A  is  located  at  the  origin  and 
sensor  B  at  (u,v),  where  u*=300  nm  (nautical  miles)  and  v-=0  nm,  in  the  system  plane.  Each 
site  is  assumed  to  have  a  range  bias  of  1  nm  and  azimuth  bias  of  0.0087  radians.  The  track 
patterns  shown  represent  typical  flight  tracks  that  may  be  encountered  by  two  radars  in  a 
practical  situation.  In  track  pattern  (a),  the  target  measurements  are  distributed  on  both 
sides  of  the  site  line  joining  sensors  A  and  B,  and  the  track  has  a  slope  of  1  relative  to  the 
site  line.  Track  (b)  is  simulated  to  be  parallel  to  the  line  connecting  the  two  sites  and  the 
measurements  are  on  only  one  side  of  the  site  line.  The  distance  between  the  track  and  the 
site  line,  d,  is  set  to  be  120  nm. 

The  means  and  normalized  MSE's  for  the  EML  and  the  LS  estimates  are  computed 
under  different  noise  assumptions.  The  deviation  of  the  measurement  noise  varies  from  0 
to  2  nm.  Performance  is  evaluated  through  Monte-Carlo  method  and  each  test  is  repeated 
100  times  to  obtain  the  averaged  results.  From  Tables  1  and  2,  we  observe  that  both  the 
EML  and  the  LS  estimates  are  unbiased.  Figures  3  and  4  show  the  variations  of  the 
normalized  MSE  of  the  registration  error  estimates  via  the  measurement  noise  deviation. 
They  clearly  demonstrate  the  improved  performance  of  the  EML  estimates  over  the  LS 
estimates,  especially  for  target  pattern  (b)  where  the  EML  algorithm  maintains  a  good 
estimation  performance  even  in  the  presence  of  strong  measurement  noise  while  the  LS 
method  performs  poorly  and  fails  as  the  measurement  noise  becomes  large. 
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Figure  5  and  6  show  the  variations  of  the  normalized  MSE's  of  the  registration 
estimates  versus  the  number  of  measurements.  The  deviation  of  the  measurement  noise 
an  is  set  to  0.5  nm,  and  the  number  of  measurements  varies  from  10  to  60  with  a  step  size 
of  10.  Again,  for  each  measurement  size,  100  trials  are  repeated  to  obtain  the  average. 
Both  EML  and  LS  estimates  converge  to  the  true  parameter  values  as  the  number  of 
measurement  increases.  However,  for  small  number  of  measurements,  the  EML  algorithm 
has  a  much  smaller  MSE  than  the  LS  method. 

The  distance  between  the  measurements  form  different  sensors  before  and  after 
registration  is  another  index  for  measuring  the  quality  of  the  registration  algorithm.  In 
Table  3,  we  compute  the  distances  between  the  original  tracks  and  between  the  updated 
tracks  using  both  the  EML  and  the  LS  methods  against  the  measurement  noise  deviations. 
Apparently,  both  methods  can  successfully  reduce  the  distance  between  the  updated  tracks 
from  different  sensors  after  registration.  The  EML  method  outperforms  the  LS  method  in 
the  sense  that  it  produces  smaller  distances  between  the  updated  tracks  from  different 
sensors  for  all  levels  of  measurement  noise  deviations.  It  should  be  mentioned  that,  unlike 
the  LS  algorithm,  the  EML  method  provides  optimal  estimates  of  the  target  tracks 
simultaneously  with  the  registration  estimates. 

To  understand  the  efficiency  of  the  EML  method  in  a  practical  environment,  real- 
life  multiple  radar  data  are  collected  from  an  air  surveillance  network.  The  radar  network 
is  composed  of  seventeen  24-by-24  foot  phased  array  L-band  long  range  radars  located 
along  the  coastline  of  Canada.  Tracks  of  air  targets  arise  from  commercial  flights.  The 
operating  specification  of  the  radars  are  summarized  as  follows: 

•  operating  frequency:  1215-1400  Mhz. 

•  instrument  range:  5  -  200  nm 

•  azimuth  coverage:  120  degrees  in  12  seconds 

•  range  resolution:  300  meters 

•  azimuth  beamwidth:  2.2  degrees 

•  probability  of  detection:  0.75 

To  remove  the  effects  of  false  targets  (clutters),  target  ID's  provided  by  the 
identification  of  friend  and  foe  (IFF)  beacon  are  used  to  extract  the  true  aircraft  tracks 
from  radar  returns  for  registration  calculation.  Each  site  in  the  system  has  two  satellite 
dishes  for  communications,  and  the  communications  between  the  radars  are  via  the  Anik 
satellite.  The  radar  network  employs  the  stereographic  projection  [14,15]  to  map  the 
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elliptic  earth  onto  a  plane  to  get  the  ground  range  from  the  slant  range  of  a  target.  The 
target  azimuth  is  measured  relative  to  the  true  north  at  the  radar  location  and  is  adjusted 
so  that  it  is  relative  to  the  true  north  at  the  origin  of  the  common  co-ordinate  system.  For 
the  data  set  used  in  this  study,  the  location  of  radar  A  is  transformed  to  the  origin  and 
hence  radar  B  is  located  at  (u,v)-(272.6406,  19.3287)  nm,  as  shown  in  Figure  7. 

Table  4  shows  the  registration  estimates  obtained  using  the  two  methods.  The 
number  of  measurements  varies  from  15  to  55  in  steps  of  10.  It  is  shown  that  the  EML 
estimates  are  more  consistent  than  the  LS  estimates  when  different  number  of 
measurement  is  used,  especially  when  this  number  is  small.  The  EML  estimates  for  K*=15 
are  close  to  those  for  K=55,  while  the  LS  method  shows  a  drastic  change  in  its  estimates. 
This  implies  that  the  EML  method  is  capable  of  satisfactory  registration  with  a  relatively 
small  number  of  measurements.  For  K-15,  the  distance  between  the  original  tracks  from 
different  sensors  is  0.7223  nm.  After  registration  using  EML,  the  distance  between  the 
updated  tracks  becomes  0.0970  nm  while  that  for  the  LS  method  is  0.6635  nm. 

Another  important  measure  of  registration  algorithm  is  its  generalization  capability. 
More  precisely,  the  registration  estimates  should  also  work  for  those  measurement  data 
not  used  in  estimating  the  registration  errors.  To  evaluate  this  aspect  of  the  two  methods, 
we  use  the  first  15  measurements  to  compute  the  system  biases.  These  estimates  are  then 
used  in  the  rest  of  the  data  set  to  see  whether  the  distance  between  the  updated  tracks  can 
be  reduced.  Figure  8  shows  the  distances  between  the  updated  tracks  versus  the  number 
of  measurements.  For  the  EML  method,  the  distances  remain  at  the  same  level  as  that  for 
the  first  15  measurements.  However,  when  the  LS  estimates  are  used,  the  generalization 
gets  worse  as  the  number  of  measurements  increases  and  the  distances  between  the 
updated  tracks  become  larger  than  that  of  the  original  tracks. 

Conclusions 

A  robust  registration  algorithm  called  the  exact  maximum  likelihood  (EML) 
method  for  multiple  sensor  data  fusion  is  developed  in  this  report.  The  EML  algorithm 
uses  the  likelihood  function  based  on  direct  measurements.  A  two-step  optimization 
procedure  is  proposed  to  estimate  the  true  target  positions  and  the  system  biases 
iteratively.  Computer  simulation  and  real  data  analysis  show  that  the  EML  method 
outperforms  the  LS  registration  algorithm.  It  can  provide  satisfactory  registration 
accuracies  with  a  small  number  of  measurements,  and  achieve  consistent  registration  error 
estimates  under  different  target  distribution  patterns.  Although  the  trade-off  for  this 
improvement  is  a  higher  computational  cost,  it  does  not  pose  any  serious  problem  for 
registration  applications  since  registration  is  an  off-line  pre-processing  procedure  for  a 
data  fusion  system. 
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Figures  and  tables 


Table  1:  Mean  of  the  registration  estimates  for  flight  track  (a) 


EML 

LS 

Vn 

A  9a 

A  rA 

A  9b 

A  7~£ 

a  eA 

A  rA 

A  &b 

A  rjg 

2.0 

0.0105 

1.1517 

0.0068 

0.8087 

0.0090 

1.0042 

0.0079 

1.0232 

bo 

0.0070 

0.9484 

0.0100 

1.0300 

0.0078 

1.1207 

0.0089 

0.9144 

1.6 

0.0117 

1.0492 

0.0066 

0.8988 

0.0084 

0.7027 

0.0096 

1.2980 

1.4 

0.0089 

0.8861 

0.0088 

0.9643 

0.0082 

0.7965 

0.0092 

1.1004 

1.2 

0.0084 

0.8543 

0.0088 

1.0838 

0.0086 

1.0502 

0.0083 

0.9232 

1.0 

0.0075 

0.8473 

0.0101 

1.1293 

0.0083 

1.0116 

0.0091 

0.9974 

0.8 

0.0094 

1.1117 

0.0084 

0.8612 

0.0084 

1.0221 

0.0091 

0.9802 

0.6 

0.0090 

1.0727 

0.0086 

0.8840 

0.0091 

1.1469 

0.0083 

0.8362 

0.4 

0.0091 

1.0857 

0.0086 

0.9275 

0.0088 

1.0765 

0.0086 

0.9618 

0.2 

0.0090 

1.0252 

0.0085 

0.9757 

0.0088 

0.9862 

0.0086 

1.0389 

0.0 

0.0087 

1.0000 

0.0087 

1.0000 

0.0087 

1.0000 

0.0087 

1.0000 

Table  2:  Mean  of  the  registration  estimates  for  flight  track  (b) 


EML 

LS 

Cn 

A  Qa 

ArA 

A  9b 

A  tb 

a  eA 

A  rA 

A  6b 

Ars 

2.0 

0.0056 

. 

1.5598 

0.0121 

1.4078 

0.0023 

2.4367 

0.0169 

1.8839 

1.8 

0.0017 

1.9003 

0.0146 

2.0610 

-0.0013 

2.2074 

0.0164 

2.6112 

1.6 

0.0060 

1.2538 

0.0112 

1.5121 

0.0049 

1.2035 

0.0109 

1.6879 

1.4 

0.0075 

0.9480 

0.0092 

1.1438 

0.0044 

1.1989 

0.0109 

1.6879 

1.2 

0.0086 

0.7097 

0.0073 

1.0020 

0.0053 

1.5171 

0.0119 

1.4776 

1.0 

0.0092 

0.6456 

0.0070 

0.9788 

0.0048 

1.1887 

0.0105 

1.7124 

0.8 

0.0092 

0.6873 

0.0073 

0.9727 

0.0061 

1.1769 

0.0102 

1.4752 

0.6 

0.0090 

0.7451 

0.0076 

0.9787 

0.0079 

1.4986 

0.0114 

1.0370 

0.4 

0.0089 

0.8183 

0.0078 

1.0190 

0.0081 

1.2191 

0.0098 

1.1019 

0.2 

0.0089 

0.8543 

0.0080 

1.0051 

0.0093 

1.0329 

0.0086 

0.8972 

0.0 

0.0087 

1.0000 

0.0087 

1.0000 

0.0087 

1.0000 

0.0087 

1.0000 
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Table  3:  Calculation  of  the  distance  between  the  updated  tracks 


TRACK  (a) 

TRACK  (b) 

Cn 

EML 

LS 

EML 

LS 

0.6004 

0.8629 

0.6142 

1.8 

0.5409 

0.8210 

0.4851 

1.6 

0.4923 

0.7412 

0.4505 

1.4 

0.3825 

0.6267 

0.3925 

1.2 

0.3177 

0.5091 

0.2988 

0.2508 

0.4273 

0.2401 

0.4605 

0.8 

0.1928 

0.3469 

0.1811 

0.6 

0.1348 

0.2619 

0.1354 

0.4 

0.0879 

0.1721 

0.0888 

0.2 

0.0434 

0.0912 

0.0440 

0.0000 

0.0000 

0.0000 

Table  4:  Registration  results  using  real  radar  data 


EML 

LS 

K 

A  eA 

ArA 

A  0B 

A  rB 

A0a 

ArA 

A8b 

Ars 

15 

-0.0046 

0.0594 

-0.0208 

-0.2895 

-0.0150 

62.6223 

0.0577 

-62.5132 

25 

-0.0069 

0.0504 

-0.0163 

-0.2177 

-0.0015 

6.4926 

-0.0148 

-6.6899 

35 

-0.0059 

-0.0158 

-0.0184 

-0.2116 

-0.0032 

1.8650 

-0.0196 

-2.0925 

45 

-0.0064 

-0.0360 

-0.0177 

-0.2035 

-0.0029 

2.9533 

-0.0185 

-3.1937 

55 

-0.0062 

-0.0675 

-0.0180 

-0.1993 

-0.0043 

1.5275 

-0.0184 

-1.7672 
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track  pattern  (a) 


site  A  site  B 


track  pattern  (b) 

Figure  2:  The  simulated  flight  track  patterns  for  registration. 
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(a)  Azimuth  Bias  Estimate  For  A  (b)  Range  Bias  Estimate  For  A 


0  0.5  1  1.5  2  0  0.5  1  1.5  2 

measurement  noise  deviation  (nm)  measurement  noise  deviation  (nm) 


(c)  Azimuth  Bias  Estimate  For  B  (d)  Range  Bias  Estimate  For  B 


0  0.5  1  1.5  2  0  0.5  1  1.5  2 

measurement  noise  deviation  (nm)  measurement  noise  deviation  (nm) 


Figure  3:  The  variation  of  the  normalised  MSE’s  of  the  registration  estimates  versus  the  measurement 
noise  deviation  for  track  pattern  (a).  The  solid  line  represents  the  EML  results,  and  the  line  labeled 
’o’  denotes  the  results  obtained  by  the  LS  method. 
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(a)  Azimuth  Bias  Estimate  For  A  (b)  Range  Bias  Estimate  For  A 


measurement  noise  deviation  (nm)  measurement  noise  deviation  (nm) 

(c)  Azimuth  Bias  Estimate  For  B  (d)  Range  Bias  Estimate  For  B 


0  0.5  1  1.5  2  0  0.5  1  1.5  2 

measurement  noise  deviation  (nm)  measurement  noise  deviation  (nm) 


Figure  4:  The  variation  of  the  normalised  MSE’s  of  the  registration  estimates  versus  the  measurement 
noise  deviation  for  track  pattern  (b).  The  solid  line  represents  the  EML  results,  and  the  line  labeled 
’o’  denotes  the  results  obtained  by  the  LS  method. 


P500711.PDF  [Page:  24  of  48] 


2  1 


(c)  Azimuth  Bias  Estimate  For  B  (d)  Range  Bias  Estimate  For 


Figure  5:  The  variation  of  the  normalised  MSE’s  of  the  registration  estimates  versus  the  number  of 
measurement  for  track  pattern  (a).  The  solid  line  represents  the  EML  results,  and  the  line  labeled  ’o’ 
denotes  the  LS  results. 


normalized  MSE  normalized  MSE 
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(a)  Azimuth  Bias  Estimate  For  A 


number  of  measurements 


Figure  6:  The  variation  of  the  normalised  MSE’s  of 
measurement  for  track  pattern  (b).  The  solid  line  rej 
denotes  the  LS  results. 


normalized  MSE  normalized  MSE 
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(b)  Range  Bias  Estimate  For  A 


number  of  measurements 


distance  (nm) 
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Appendix:  Matlab  Program  Listing 

rgdatal.m 

%  This  function  file  rgdataLm  generates  the  simulated  target  data.  The 
%  returns  are  azimuth  and  range  with  respect  to  site  A  and  B,  respectively. 
%  The  inputs  are  site  B  coordinate,  the  slope,  number  of  target  K,  the  step 
%  size,  the  pattern  number  and  myv  is  the  y  functin  value  for  pattern  3  and 
%  4. 

%  Yifeng  Zhou,  McMaster  University,  Jan.  1995 

function  rgm=rgdatal(u,K,st) 

rgm-zeros(4,K); 

mta=»zeros(l,K); 

mra«zeros(l,K); 

mtb“zeros(l,K); 

mrb=zeros(l,K); 

my-zeros(l,K); 

mx-u*ones(  1  ,K)/2; 
mk=0:K-l; 
my-(mk-(K- 1  )/2)*  st; 


thrd-10e-4; 


for  k-l:K 

mra(k)=sqrt(mx(k)*mx(k)+my(k)*my(k»; 

mrb(k)=sqrt((mx(k)-u)*(mx(k)-u)+my(k)*my(k)); 

if  abs(my(k))<thrd 
mta(k)=pi/2; 
mtb(k)=-pi/2; 
else 

if  my(k)  >=  0 

mta(k)=atan(mx(k)/my(k)); 

mtb(k)=atan((mx(k)-u)/my(k)); 

else 

mta(k)=-atan(my(k)/mx(k))+pi/2; 

mtb(k)=-atan(my(k)/(mx(k)-u))-pi/2; 

end 

end 

aid 


rgm“[mta;mra;mtb;mrb] ; 
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rgdata2.m 

%  This  function  file  rgdata2.m  generates  the  simulated  target  data  for 
%  pattern  2.  The  returns  are  azimuth  and  range  with  respect  to  site  A 
%  and  B,  respectively.  The  inputs  are  site  B  coordinate,  the  slope, 

%  number  of  target  K,  the  step  size. 

%  Yifeng  Zhou,  McMaster  University,  Jan.  1995 

function  rgm=rgdata2(u,kr,K,st) 

rgm=zeros(4,K); 

mta=zeros(l,K); 

mra=zeros(l,K); 

mtb=zeros(l,K); 

mrb=zeros(l,K); 

mx-zeros(l,K); 

my=zeros(l,K); 


mk=l:K; 

mx=(mk- 1  -(K- 1  )/2)  *st+u/2; 
my=kr*(mx-u/2); 

thrd—10e-4; 


for  k«=l:K 

mra(k)=sqrt(mx(k)  *mx(k)+my  Ok.)  *my(k)); 
mrb(k)=sqrt((mx(k)-u)*(mx(k)-u)+my(k)*my(k)); 

if  abs(my(k))<thrd 
mta(k)=pi/2; 
mtb(k)=-pi/2; 
else 

if  my(k)  >=  0 

mta(k)=atan(mx(k)/my(k)); 
mtb(k)=atan((mx(k)-u)/  my(k)) ; 
else 

mta(k)=-atan(my(k)/mx(k))+pi/2; 

mtb(k)=-atan(my(k)/(mx(k)-u))-pi/2; 

end 

end 

end 


rgm*-[mta;mra;mtb;mrb]; 
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rgdata3.ni 

%  This  function  file  rgdata3.m  generates  the  simulated  target  data  for 
%  pattern  3.  The  returns  are  azimuth  and  range  with  respect  to  site  A 
%  and  B,  respectively.  The  inputs  are  site  B  coordinate,  the  slope, 

%  number  of  target  K,  step  size  and  myv  is  the  y  functin  value. 

%  Yifeng  Zhou,  McMaster  University,  J an.  1 995 

function  rgm=rgdata3(u,K,st,myv) 


rgm=zeros(4,K); 

mta=*zeros(l,K); 

mra-zeros(l,K); 

mtb=zeros(l,K); 

mrb=zeros(l,K); 

mx=zeros(l,K); 

my«=zeros(l,K); 


mk-l:K; 

mx=(mk- 1  -(K- 1  )/2)*st+u/2; 
my=-myv*ones(l  ,K); 
thrd-lOe-4; 


fork=l:K 

mra(k)=sqrt(mx(k)*mx(k)+my(k)*my(k)); 

mrb(k)=sqrt((mx(k)-u)*(mx(k)-u)+my(k)*my(k)); 

if  abs(my(k))<thrd 
mta(k)=pi/2; 
mtb(k)*=-pi/2; 
else 

if  my(k)  >-  0 

mta(k)=atan(mx(k)/my(k)); 

mtb(k)=atan((mx(k)-u)/my(k»; 

else 

mta(k)=-atan(my(k)/  mx(k))+pi/2; 
mtb(k)=-atan(my(k)/(mx(k)-u))-pi/2; 

end 

end 

end 


rgm=[mta;mra;mtb;mrb]; 
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rgdata4.m 

%  This  function  file  rgdata4.m  generates  the  simulated  target  data  for 
%  pattern  4.  The  returns  are  azimuth  and  range  with  respect  to  site  A 
%  and  B,  respectively.  %  The  inputs  are  site  B  coordinate,  number  of 
%  target  K,  the  step  size  and  the  y  function  value  myv. 

%  Yifeng  Zhou,  McMaster  University,  Jan.  1995 

function  rgm=rgdata4(u,K,st,myv) 

rgm=zeros(4,K); 

mta=zeros(l,K); 

mra=zeros(l,K); 

mtb=zeros(l,K); 

mrb=zeros(l,K); 

mx=zeros(l,K); 

my-zeros(l,K); 

HK=K/2; 

mk=l:HK; 

mx(  1 :  HK)=(mk- 1  -(HK- 1  )/2)  *s  t+u/2; 
mx(HK+ 1 :  K)=mx(  1 :  HK); 
my(l  :HK)=myv*ones(  1  ,HK); 
my(HK+l:K)=-myv*ones(l,K-HK); 
thrd=10e-4; 

for  k-l:K 

mra(k)=sqrt(mx(k)*mxOc)+my(k)*my(k)); 

mrb(k)=sqrt((mx(k)-u)*(mx(k)-u)+my(k)*my(k)); 

if  abs(my(k))<thrd 
mta(k)=pi/2; 
mtb(k)=-pi/2; 
else 

if  my(k)  >=  0 

mta(k)=*atan(mx(k)/my  (k)) ; 
mtb(k)=atan((mx(k)-u)/my  (k)); 
else 

mta(k)=-atan(my(k)/mx(k))+pi/2; 
mtb(k)=-atan(my(k)/  (mx(k)-u))-pi/2; 
end 
end 
end 


rgm= [mta;  mra;  mtb;  mrb] ; 
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rgsaml.m 

%  This  file  rgsamLm  is  is  to  compute  the  MSE  of  system  error  via  the  changes 
%  of  number  of  measurements  for  the  EML  and  the  LS  approaches. 

cf0  Yifeng  Zhou,  McMaster  University  1995 


clear; 

kn- 10: 10:60; 
lkn=length(kn); 

%  dta,  dra,  dtb,drb  denotes  the  system  bias  of  site  A  and  B 

dta=0.5*pi/180; 

dtb=0.5*pi/180; 

dra=1.0; 

drb-1.0; 

etaO=[dta,dra,dtb,drb]'; 

%  the  unit  of  range  measurement  is  nautical  miles 

u=300; 

v=0; 


cmeta=zeros(4,lkn); 

lmeta~zeros(4,lkn); 

emeta=zeros(4,lkn); 

laeta-=zeros(4,lkn); 

eaeta=»zeros(4,lkn); 

eiter-=zeros  ( 1  ,lkn) ; 

ldse=zeros(l  ,lkn); 
edse=zeros(l  ,lkn); 

TM=100; 


thrd=1.0e-3; 

mu0=0.8; 

epsilone=0.02; 

epsilonx=0.01; 

tumax-21; 

tu=»zeros(l  ,tumax); 

tmpxi=zeros(2,tumax); 

mu“[mu0.A(l:tumax-l),0]; 
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st«=4; 

mgv-120; 

kr=l; 

snr=0.5*dra; 

t^***************  compute  the  asymptotic  covariance  matrix  ******* 

Rk=zeros(4,2); 

tceta=*zeros(4,4); 

ceta=zeros(4,4); 

for  tn=l:lkn 

K=kn(tn); 

mta=zeros(l,K); 

mtb=«zeros(l,K); 

mra=zeros(l,K); 

mrb=zeros(l,K); 

tm=rgdatal  (u,K,st); 

%tm-rgdata2(u,kr,K,st); 

%tm=rgdata3(u,K,st,mgv); 

%tm=rgdata4(u,K,st,mgv); 

mta=tm(l,:); 

mra=tm(2,:); 

mtb=tm(3,:); 

mrb=tm(4,:); 


%  simulate  the  biased  xy  coordinate  with  respect  to  site  A  and  B 

txa-=zeros(l,K); 

tya=zeros(l,K); 

txb=zeros(l,K); 

tyb=zeros(l,K); 

txa=(mra+dra).*sin(mta+dta); 
tya-=(mra+dra).*cos(mta+dta); 
txb=(mrb+drb).*sin(mtb+dtb); 
tyb=(mrb+drb).*cos(mtb+dtb) ; 

for  k=l:K 

alpha=mra(k)*  sin(mta(k)); 
beta=mra(k)*cos(mta(k)); 
rda=sqrt(alphaA2+betaA2); 
rdb=sqrt((alpha-u)A2+(beta-v)A2); 
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Rk(  1,1)— 1  -eta0(2)  *betaA2/rdaA3; 

Rk(2,l)-=etaO(l)+etaO(2)*alpha*beta/rdaA3; 

Rk(3,l)—l-etaO(4)*(beta-v)A2/rdbA3; 

Rk(4, 1 ) =eta0(3  )+etaO  (4)  *  (alpha-u)  *  (beta-v)/rdb  A3 ; 

Rk(l,2)=-etaO(l)+etaO(2)*alpha*beta/rdaA3; 

Rk(2,2)=- 1  -etaO(2)*alphaA2/rdaA3; 
Rk(3,2)=-etaO(3)+etaO(4)*(alpha-u)*(beta-v)/rdbA3; 
Rk(4,2)=- 1  -etaO(4)*(alpha-u)A2/rdbA3 ; 

Ppk-eye^^^^invCRk^Rk^Rk'; 

AA=zeros(4,4); 

AA(1 :2, 1 :2)=[beta,alpha/rda;-alpha,beta/rda]; 
AA(3:4,3:4)=[beta-v,(alpha-u)/rdb;-(alpha-u),(beta-v)/rdb]; 

tceta=tceta+AA'*Ppk*AA; 

end 

ceta=inv(tceta); 

cmeta(:  ,tn)=snr*snr*diag(ceta); 


%  *********  obtain  the  initial  estimate  on  xi  and  eta  *********** 

txit«*zeros(2,K); 

xi=zeros(2,K); 

xi(l,:)=mra.*sin(mta); 

xi(2,:)=mra.*cos(mta); 

xiO=xi; 

tda=sqrt(xi(l,:).A2+xi(2,:).A2); 

tdb=sqrt((xi(l,:)-u).A2+(xi(2,:)-v).A2); 

for  tms"l:TM 

%  generate  the  simulated  target  data  for  different  SNR 

xa=txa+snr  *randn(  1  ,K); 

ya=tya+snr*randn(l  ,K); 

xb=txb+snr*randn(  1  ,K)+u; 

yb=tyb+snr*randn(l  ,K)+v; 

X-[xa;ya;xb;yb]; 

XA-[xa;ya]; 

XB-[xb;yb]; 

%  obtain  the  noise  azimuth  and  range  measurement 
for  k-l:K 

ra(k)=sqrt(xa(k)  *xa(k)+ya(k)*ya(k)); 
rb(k)=sqrt((xb(k)-u)*(xb(k)-u)+(yb(k)-v)*(yb(k)-v»; 
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if  abs(ya(k))<thrd 
ta(k)=pi/2; 
elseif  ya(k)>=0 
ta(k)=atan(xa(k)/ya(k)); 
else 

ta(k)-=-atan(yaGc)/xa(k))+pi/2; 

end 

if  abs(yb(k))<thrd 
tb(k)=-pi/2; 
elseif  yb(k)>=0 
tb(k)=atan((xb(k)-u)/yb(k)); 
else 

tb(k)=-atan(yb(k)/(xb(k)-u))-pi/2; 

end 

end 

%  hc***************  the  least  squares  algorithm  ******************** 

LA-zeros(2*K,4); 

bl=zeros(2*K,l); 

for  k=l:K 
ml-2*(k-l)+l; 
m2-2*k; 

LA(ml,:)=[ra(k)*cos(ta(k)),sin(ta(k)),-rb(k)*cos(tb(k)),-sin(tb(k))]; 

LA(m2,:)=[-ra(k)*sin(ta(k)),cos(ta(k)),rb(k)*sm(tb(k)),-cos(tb(k))]; 

bl(ml ,  l)=xa(k)-xb(k); 
bl(m2, 1  )=*ya(k)-yb(k); 
end 

leta=inv(LA'*LA)*LA'*bl; 

txit(  1  ,:)=(ra-leta(2)).*sin(ta-leta(  1 )); 

txit(2,:  )=(ra-leta(2» .  *cos(ta-leta(  1 )); 

tdsel-=(txit(l,:)-xiO(l  ,:))*(txit(l  ,:)-xiO(l 

tdse2=(txit(2,:)-xi0(2,:))*(txit(2,:)-xi0(2,:))'; 

tdse=sqrt(tdse  1 +tdse2)/K; 

ldse(tn)=ldse(tn)+tdse/TM; 

%  **********  the  exact  maximum  likelihood  algorithm  ****************** 
TAl=zeros(2,2); 

TA2=*zeros(2,2); 

tal=zeros(2,l); 

ta2=zeros(2,l); 
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%  evaluation  of  the  system  error  eta 
for  k=l:K 
alpha=xi(l,k); 
beta=xi(2,k); 
rda-tda(k); 
rdb=tdb(k); 

A 1 = [beta,alpha/rda;  -alpha, beta/rda] ; 
A2-[beta-v,(alpha-u)/rdb;-(alpha-u),(beta-v)/rdb]; 

TA1-TA1+A1'*A1; 
tal-tal+Al '*(XA(:  ,k)-xi(:,k)); 

TA2-TA2+A2'*A2; 
ta2«=ta2+A2'*(XB(:  ,k)-xi(:  ,k)); 
end 

eta=[inv(TAl)*tal;inv(TA2)*ta2]; 

%  estimate  the  true  target  position  xi,  Gr  is  the  gradient 
Dm=zeros(4,l); 

Gr-zeros(2,l); 

Ja»zeros(2,4); 

Hes=zeros(2,2); 

for  k-l:K 
alpha-xi(l,k); 
beta-xi(2,k); 
rda=tda(k); 
rdb=tdb(k); 

Ja(  1 , 1  )=- 1  -eta(2)  *betaA2/rdaA3 ; 

Ja(  1 ,2)=eta(l  )+eta(2)*alpha*beta/rdaA3; 

Ja(l,3)~  l-eta(4)*(beta-v)A2/rdbA3; 

Ja(l,4)=eta(3)+eta(4)*(alpha-u)*(beta-v)/rdbA3; 

Ja(2,l)=-eta(l)+eta(2)*alpha*beta/rdaA3; 

Ja(2,2)— l-eta(2)*alphaA2/rdaA3; 
Ja(2,3)=-eta(3)+eta(4)*(alpha-u)*(beta-v)/rdbA3; 
Ja(2,4)— l-eta(4)*(alpha-u)A2/rdbA3; 

Dm(  1  )=alpha+alpha*eta(2)/rda+beta*eta(  1 ); 
Dm(2)=beta+beta*eta(2)/rda-alpha*eta(l); 
Dm(3)=alpha+(alpha-u)*eta(4)/rdb+(beta-v)*eta(3); 
Dm(4)=beta+(beta-v)*eta(4)/rdb-(alpha-u)*eta(3); 

Gr-2*  Ja*(X(:  ,k)-Dm); 

Hes=2*Ja*Ja'; 

Mdn=inv(Hes)*Gr; 
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hmdn-=sqrt(Mdn'*Mdn); 

while  hmdn>epsilonx 

%  determine  the  step  in  the  gradient 
for  mx=l:tumax 

tmpxi(:  ,xnx)=xi(:  ,k)-mu(mx)*Mdn; 
alpha=tmpxi(  1  ,mx); 
beta=tmpxi(2,mx); 
rda-=sqrt(alphaA2+betaA2); 
rdb=sqrt((alpha-u)A2+(beta-v)A2); 

Dm(  l)=alpha+alpha*eta(2)/rda+beta*eta(l ); 
Dm(2)=beta+beta*eta(2)/rda-alpha*eta(l); 
Dm(3)=alpha+(alpha-u)*eta(4)/rdb+(beta-v)*eta(3); 
Dm(4)=beta+(beta-v)*eta(4)/rdb-(alpha-u)’i<eta(3); 


tu(mx)=(X(:  ,k)-Dm)'*(X(:  ,k)-Dm); 
end 

[tumx,tumd]=min(tu); 
xi(:  ,k)=tmpxi(:  ,tumd); 
alpha=xi(l,k); 
beta=xi(2,k); 

rda=sqrt(alphaA2+betaA2); 

rdb=sqrt((alpha-u)A2+(beta-v)A2); 

tda(k)=rda; 

tdb(k)=rdb; 

Ja(  1,1  )=- 1  -eta(2)  *betaA2/rdaA3; 

Ja(l,2)=eta(l)+eta(2)*alpha*beta/rdaA3; 

Ja(l,3)=-l-eta(4)*(beta-v)A2/rdbA3; 

Ja(l,4)=eta(3)+eta(4)*(alpha-u)*(beta-v)/rdbA3; 

Ja(2,l)=-eta(l)+eta(2)*alpha*beta/rdaA3; 

Ja(2,2)*=- 1  -eta(2)*alphaA2/rdaA3; 

Ja(2,3)=-eta(3)+eta(4)*(alpha-u)*(beta-v)/rdbA3; 

Ja(2,4)=-l-eta(4)*(alpha-u)A2/rdbA3; 

Dm(  1  )=alpha+alpha*eta(2)/rda+beta*eta(l ); 
Dm(2)=beta+beta*eta(2)/rda-alpha*eta(l); 
Dm(3)=alpha+(alpha-u)*eta(4)/rdb+(beta-v)*eta(3); 
Dm(4)=beta+(beta-v)*eta(4)/rdb-(alpha-u)*eta(3); 

Gr==2*Ja*(X(:,k)-Dm); 

Hes»=2*Ja*Ja'; 

Mdn=inv(Hes)*Gr; 
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hmdn=sqrt(Mdn'  *Mdn); 
end 
end 

TAl-zeros(2,2); 

TA2-zeros(2,2); 

tal=zeros(2,l); 

ta2«=zeros(2,l); 

for  k-=l:K 
alpha“xi(l,k); 
beta=xi(2,k); 
rda-=tda(k); 
rdb=tdb(k); 

Al=[beta,alpha/rda; -alpha, beta/rda] ; 
A2-[beta,(alpha-u)/rdb;-(alpha-u),(beta-v)/rdb]; 

TA1-TA1+A1'*A1; 
ta  1 =ta  1  +  A 1 '  *  (XA(:  ,k)-xi(:  ,k)); 
TA2-TA2+A2'*A2; 
ta2*=ta2+A2'*(XB(:  ,k)-xi(:  ,k)); 
end 

teta-[inv(TAl)*tal;inv(TA2)*ta2]; 
adis«sqrt((eta-teta)'*  (eta-teta)) ; 
eta-teta; 

emiter=l; 

%  begin  the  iteration  process 
while  adis>epsilone 

emiter=emiter+l; 

for  k— 1:K 
alpha=xi(l,k); 
beta-xi(2,k); 
rda=tda(k); 
rdb=tdb(k); 

Ja(  1 , 1  )=- 1  -eta(2)  *betaA2/  rdaA3 ; 

Ja(  1 ,2)=eta(  1  )+eta(2)*alpha*beta/rdaA3; 
Ja(l,3)-=-l-eta(4)*(beta-v)A2/rdbA3; 
Ja(l,4)=eta(3)+eta(4)*(alpha-u)*(beta-v)/rdbA3; 
Ja(2, 1  )=-eta(  1  )+eta(2)*alpha*beta/rdaA3; 
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J  a(2,2)=- 1  -eta(2)  *alphaA2/  rdaA3; 

Ja(2,3)=-eta(3)+eta(4)*(alpha-u)*(beta-v)/rdbA3; 
Ja(2,4)=- 1  -eta(4)*(alpha-u)A2/rdbA3; 

Dm(  1  )=alpha+alpha*  eta(2)/rda+beta*eta(l ); 
Dm(2)=beta+beta*eta(2)/  rda-alpha*eta(  1 ); 

Dm(3) =alpha+ (alpha-u)  *  eta(4)/  rdb+(beta-v )  *  eta(3); 

Dm(4)“beta+(beta-v)*eta(4)/rdb-(alpha-u)*eta(3); 

Gr=2*Ja*(X(:,k)-Dm); 

Hes=2*Ja*Ja'; 

Mdn=-inv(Hes)*Gr; 
hmdn=sqrt(Mdn'  *Mdn); 

while  hmdn>epsilonx 
%  determine  the  step  in  the  gradient 
for  mx=l:tumax 

tmpxi(:  ,mx)=xi(:  ,k)-mu(mx)>|cMdn; 
alpha=tmpxi(  1  ,mx); 
beta=tmpxi(2,mx); 
rda=sqrt(alphaA2+betaA2); 
rdb=sqrt((alpha-u)A2+(beta-v)A2); 

Dm(l)=alpha+alpha*eta(2)/rda+beta*eta(l); 
Dm(2)=beta+beta*eta(2)/rda-alpha*eta(  1 ); 
Dm(3)=alpha+(alpha-u)*eta(4)/rdb+(beta-v)*eta(3); 
Dm(4)=*beta+(beta-v)*eta(4)/rdb-(alpha-u)*eta(3); 


tu(mx)=(X(:,k)-Dm)'*(X(:,k)-Dm); 

end 

[tumx,tumd]=min(tu); 
xi(:  ,k)=tmpxi(:  ,tumd); 
alpha=xi(l,k); 
beta=xi(2,k); 

rda=sqrt(alphaA2+betaA2); 

rdb=sqrt((alpha-u)A2+(beta-v)A2); 

tda(k)=rda; 

tdb(k)=rdb; 

Ja(  1 , 1  )=- 1  -eta(2)*betaA2/rdaA3; 

Ja(  1 ,2)=eta(  1  )+eta(2)*alpha*beta/rdaA3; 
Ja(l,3)=-l-eta(4)*(beta-v)A2/rdbA3; 

Ja(  1 ,4)=eta(3)+eta(4)*(alpha-u)*(beta-v)/rdbA3; 
Ja(2, 1  )=-eta(  1  )+eta(2)*alpha*beta/rdaA3; 

J a(2,2)=- 1  -eta(2)  *  alpha A2/rdaA3 ; 
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J  a(2,3)=-eta(3)+eta(4)  *  (alpha-u)  *  (beta-v )/ rdbA3 ; 
Ja(2,4)=- 1  -eta(4)  *  (alpha-u)  A2/rdbA3; 

Dm(  1  )-*alpha+alpha*eta(2)/  rda+beta*  eta(  1  )*, 
Dm(2)=beta+beta*eta(2)/rda-alpha*eta(l); 
Dm(3)-=alpha+(alpha-u)*eta(4)/rdb+(beta-v)*eta(3); 
Dm(4)=beta+(beta-v)*eta(4)/rdb-(alpha-u)*eta(3); 

Gr=2* J  a*  (X(:  ,k)-Dm); 

Hes«2*Ja*Ja'; 

Mdn*=inv(Hes)*Gr, 

hmdn=sqrt(Mdn'*Mdn); 

end 

end 

TAl«zeros(2,2); 

TA2=zeros(2,2); 

tal«zeros(2,l); 

ta2»zeros(2,l); 

for  k-l:K 
alpha*»xi(l,k); 
beta*=xi(2,k); 
rda=tda(k); 
rdb=tdb(k); 

A 1  ■=  [beta, alpha/  rda;-alpha,beta/  rda] ; 

A2»[beta,(alpha-u)/rdb;-(alpha-u),(beta-v)/rdb]; 

TA1-TA1+A1'*A1; 

tal=tal+Al'*(XA(:,k)-xi(:,k)); 

T  A2“T  A2+A2'  *  A2; 
ta2»=ta2+A2,*(XB(:  ,k)-xi(:  ,k)); 

end 

teta=[inv(TAl)*tal;inv(TA2)*ta2]; 

adis=sqrt((eta-teta)'*(eta-teta)); 

eta=teta; 

end 

eiter(tn)=eiter(tn)+emiter/TM; 

lmeta(:  ,tn)=lmeta(:  ,tn)+(leta-etaO)  .A2/TM; 
emeta(:  ,tn)=emeta(:  ,tn)+(eta-etaO).A2/TM; 
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laeta(:  ,tn)=laeta(:  ,tn)+leta/TM; 
eaeta(:  ,tn)=eaeta(:  ,tn)+eta/TM; 

tdse  l»=(xi(  1  ,:)-xi0(  1 ))*(xi(l  ,:)-xi0(l 
tdse2=(xi(2,:)-xi0(2,:))*(xi(2,:)-xi0(2,:))'; 
tdse=sqrt(tdse  1 +tdse2)/K; 
edse(tn)=edse(tn)+tdse/TM; 

end 


save  rgsaml.mat; 
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drgsnrl.m 

%  This  file  drgsnrl.m  is  is  to  compute  the  MSE  and  mean  of  the  system  biases 
%  via  the  changes  of  SNR  for  1)  the  exact  maximum  likelihood  registration 
%  algorithm  2)  the  least  squars  algorithm 

%  McMaster  University  1995 

clear, 

global  u  v  geta  globalx; 

%  K  is  the  number  of  targets  used  in  the  registration  algorithm 
K=20; 

%  dta,  dra,  dtb.drb  denotes  the  system  bias  of  site  A  and  B 

dta-0.5*pi/180; 

dtb=0.5*pi/180; 

dra=l; 

drb=l; 

etaO=[dta,dra,dtb,drb]'; 

%  the  unit  of  range  measurement  is  nautical  miles 

u=300; 

v=0; 


%  obtain  the  true  azimuth  and  range  information  with  respect  to  A  and  B 

%  kr  is  the  slope,  st  is  the  step  size,  mgv  is  the  function  for  pattern 

%  "3"  and  "4",  pd  denotes  the  pattern  number 

mta=zeros(l,K); 

mtb=zeros(l,K); 

mra=zeros(l,K); 

mrb=zeros(l,K); 


st=4; 

mgv- 120; 
kr-1; 


tm=rgdata  1  (u,K,st); 
%tm=rgdata2(u,kr,K,st); 
%tm=*rgdata3(u,K,st,mgv); 
%tm-rgdata4(u,K,st,mgv); 


mta=tm(l,:); 

mra-tm(2,:); 

mtb=tm(3,:); 

mrb=tm(4,:): 
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%  simulate  the  biased  xy  coordinate  with  respect  to  site  A  and  B 

txa=zeros(l,K); 

tya=zeros(l,K); 

txb=zeros(l,K); 

tyb=zeros(l,K); 


txa*=(mra+dra) .  *sin(mta+dta); 

tya=(mra+dra).*cos(mta+dta); 
txb=(mrb+drb)  .*  sin(mtb+dtb); 
tyb=(mrb+drb).*cos(mtb+dtb); 


snr0=0*dra; 

snrl=2.0*dra; 

dsnr=0.2*dra; 

snr~snrl:-dsnr:snrO; 

lsnr=-length(snr); 

cmeta=zeros(4,lsnr); 

lmeta=*zeros(4,lsnr); 

emeta=zeros(4,lsnr); 

dmeta«zeros(  1  ,lsnr) ; 
dleta*=zeros(l  ,lsnr); 

aeeta-zeros  (4,lsnr); 
aleta=zeros(4,lsnr); 

eiter=zeros(l  ,lsnr); 

thrd=1.0e-4; 
epsilone=0.001 ; 

%***************  compute  the  asymptotic  covariance  matrix 

Rk=zeros(4,2); 

tceta=zeros(4,4); 

ceta=zeros(4,4); 

for  k*=l:K 

alpha=mra(k)*sin(mta(k)); 

beta=mra(k)*cos(mta(k)); 

rda=sqrt(alphaA2+betaA2); 

rdb=sqrt((alpha-u)A2+(beta-v)A2); 


******** 
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Rk(l  ,1)«=-  l-etaO(2)*betaA2/rdaA3; 

Rk(2,  l)»etaO(  l)+etaO(2)*alpha*beta/rdaA3; 

Rk(3, 1)=-  l-etaO(4)*(beta-v)A2/rdbA3; 
Rk(4,l)=etaO(3)+etaO(4)*(alpha-u)*(beta-v)/rdbA3; 

Rk(l  ,2)=-eta0(l  )+etaO(2)*alpha*beta/rdaA3; 

Rk(2,2) =- 1  -etaO  (2)  *  alphaA2/rdaA3 ; 
Rk(3,2)=-etaO(3)+etaO(4)*(alpha-u)*(beta-v)/rdbA3; 
Rk(4,2)—  1  -eta0(4)  *(alpha-u)  A2/rdbA3 ; 

Ppk^eye(4,4)-Rk*inv(Rk'*Rk)*Rk'; 

AA-zeros(4,4); 

AA(l:2,l:2)=»[beta,alpha/rda;-alpha,beta/rda]; 

AA(3 : 4,3 : 4)= [b  eta-v ,  (alpha-u)/ rdb;  -  (alpha-u)  ,(beta-v  )/rdb] ; 

tceta=tceta+ AA'  *Ppk*  AA; 
end 

ceta=inv(tceta); 
for  tn=l:lsnr 

cmeta(:  ,tn)=snr(tn)*snr(tn)*diag(ceta); 
end 


%  *********  obtain  the  initial  estimate  on  xi  and  eta  *********** 

xiO(  1 )=mra.  *  s  in(mta) ; 

xi0(2,:)=mra.*cos(mta); 

TM=100; 

for  tn=l:lsnr 

for  tms=l:TM 

%  generate  the  simulated  target  data  for  different  SNR 
xa-=txa+snr(tn)*randn(l  ,K); 
ya=tya+snr(tn)*randn(l  ,K); 
xb=txb+snr(tn)*randn(l  ,K)+u; 
y  b=tyb+snr  (tn)  *randn(  1  ,K)+v; 

X=[xa;ya;xb;yb]; 

XA«[xa;ya]; 

XB=[xb;yb]; 


%  obtain  the  noise  azimuth  and  range  measurement 
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for  k=»l:K 

ra(k)=sqrt(xa(k)*xa(k)+ya(k)*ya(k)); 

rb(k)=sqrt((xb(k)-u)*(xb(k)-u)+(yb(k)-v)*(yb(k)-v)); 

if  abs(ya(k))<thrd 
ta(k)=pi/2; 
tb(k)=-pi/2; 
elseif  ya(k)>=0 
ta(k)-atan(xa(k)/ ya(k)); 
tb(k)=atan((xb(k)-u)/yb(k)); 
else 

ta(k)=-atan(ya(k)/xa(k))+pi/2; 
tb(k)=-atan(yb(k)/  (xb(k)-u))-pi/2; 
end 
end 

%  ****************  thg  least  squares  algorithm  ******************** 

LA-*zeros(2*K,4); 

bl=zeros(2*K,l); 

for  k«l:K 
ml-2*(k-l)+l; 
m2=2*k; 

LA(ml,:)=[ra(k)*cos(ta(k)),sin(taOc)),-rb(k)*cos(tb(k)),-sin(tbCk))]; 
LA(m2,:  )=[-ra(k)  *  sin(ta(k)),cos(ta(k)),rb(k)*  sin(tb(k)),-cos  (tb(k))] ; 
bl(m  1 , 1  )=xa(k)-xb(k); 
bl(m2,l  )=ya(k)-yb(k); 
end 

leta-inv(LA,*LA)*LA'*bl; 


%  **********  the  exact  maximum  likelihood  algorithm  ****************** 
eta=etaO; 

%fork-l:K 

%  tmpax=(ra(k)-eta(2))*sin(ta(k)-eta(l)); 

%  tmpay«(ra(k)-eta(2))*cos(ta(k)-eta(l )); 

%  tmpbx”=(rb(k)-eta(4))*sin(tb(k)-eta(3))+u; 

%  tmpby*=(rb(k)-eta(4))*cos(tb(k)-eta(3))+v; 

%  xi(l,k)=(tmpax+tmpbx)/2; 

%  xi(2,k)=(tmpay+tmpby)/2; 

%end 

for  k=l:K 
geta=eta; 
globalx=X(:,k); 
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ta2«zeros(2,l); 
for  k=l:K 

alpha-xi(l.k); 

beta=xi(2,k); 

rda-sqrt(aiphaA2+betaA2); 

rdb=sqrt((alpha-u)A2+(beta-v)A2); 

A1 =[beta,alpha/rda;-alpha,beta/rda] ; 
A2-[beta-v,(alpha-u)/rdb;-(alpha-u),(beta-v)/rdb]; 

TA1-TA1+A1'*A1; 
tal»tal+Al'*(XA(:  ,k)-xi(:  ,k)); 
TA2-TA2+A2'*A2; 
ta2=ta2+A2'*(XB(:  ,k)-xi(:  ,k)); 
end 

teta=[inv(TAl)*tal;inv(TA2)*ta2]; 

adis-sqrt((eta-teta)'*(eta-teta)); 

eta~teta; 

end 

lmeta(:,tn)=lmeta(:,tn)+Geta-etaO)  A2/TM; 
emeta(:,tn)=emeta(:,tn)+(eta-etaO)  A2/TM; 

dtxKxi(l,:)-xiO(l,:))*(xi(l,:)-xiO(l,:))'; 

dty-(xi(2,:)-xi0(2,:))*(xi(2,:)-xi0(2,:))'; 


ltxa-(ra-leta(2)) .  *  s  in(ta-leta(  1 )) ; 
ltya«=(ra-leta(2)).*cos(ta-leta(l)); 
ltxb=(rb-leta(4)).*sin(tb-leta(3))+u; 
ltyb=(rb-leta(4)).*cos(tb-leta(3))+v; 

dlxa=(ltxa-xiO(  1 , : ))  *(ltxa-xiO(  1 
dlya=(ltya-xi0(2,:))*(ltya-xi0(2,:))'; 
dlxb=(Itxb-xiO(l  ,:))*(ltxb-xiO(l 
dlyb=(ltyb-xiO(2,:))*(ltyb-xiO(2,:»'; 

dleta(tn)=dleta(tn)+sqrt((dlxa+dlxb+dlya+dlyb)/(2*TM)); 

dmeta(tn)=dmeta(tn)+sqrt((dtx+dty)/TM); 

aleta(:  ,tn)=aleta(:  ,tn)+leta/TM; 
aeeta(:  ,tn)=aeeta(:  ,tn)+eta/TM; 
end 

save  drgsnrl.mat; 
end 
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